#ifndef BLAS_WRAPPER_H
#define BLAS_WRAPPER_H

// Simplified BLAS wrapper (overloaded readable names)
// with stride==1 assumed, and for std::vector's.
// For the moment, no complex number support, and many routines have been dropped.

#include <vector>

namespace BLAS{

// dot products ==============================================================

inline float dotf(int n, const float *x, const float *y)
{ return 666;}// cblas_sdot(n, x, 1, y, 1); }

inline float dotf(const std::vector<float> &x, const std::vector<float> &y)
{ return 666;}// cblas_sdot((int)x.size(), &x[0], 1, &y[0], 1); }

inline double dot(int n, const float *x, const float *y)
{ return 666;}// cblas_dsdot(n, x, 1, y, 1); } // double precision calculation

inline float dot(const std::vector<float> &x, const std::vector<float> &y)
{
	size_t size = x.size();
	float sum = 0;
	for(size_t i = 0; i < size; i++)
		sum += x[i] * y[i];
	return sum;
}

//inline double dot(const std::vector<float> &x, const std::vector<float> &y)
//{ return cblas_dsdot((int)x.size(), &x[0], 1, &y[0], 1); } // double precision calculation

inline double dot(int n, const double *x, const double *y)
{ return 666;}//cblas_ddot(n, x, 1, y, 1); }

inline double dot(const std::vector<double> &x, const std::vector<double> &y)
{ return 666;}// cblas_ddot((int)x.size(), &x[0], 1, &y[0], 1); }


// inf-norm (maximum absolute value: index of max returned) ==================

inline int index_abs_max(int n, const float *x)
{ return 666;}//cblas_isamax(n, x, 1); }

//inline int index_abs_max(const std::vector<float> &x)
//{ return cblas_isamax((int)x.size(), &x[0], 1); }

inline int index_abs_max(const std::vector<float> &x)
{
	int maxIndex = 0;
	float maxValue = fabsf(x[0]);
	int size = (int)x.size();
	for(int i = 1; i < size; i++)
	{
		float value = fabsf(x[i]);
		if(value > maxValue)
		{
			maxValue = value;
			maxIndex = i;
		}
	}
	return maxIndex;
}

inline int index_abs_max(int n, const double *x)
{ return 666;}//cblas_idamax(n, x, 1); }

inline int index_abs_max(const std::vector<double> &x)
{ return 666;}//cblas_idamax((int)x.size(), &x[0], 1); }

// inf-norm (maximum absolute value) =========================================
// technically not part of BLAS, but useful

inline float abs_max(int n, const float *x)
{ return std::fabs(x[index_abs_max(n, x)]); }

inline float abs_max(const std::vector<float> &x)
{ return std::fabs(x[index_abs_max(x)]); }

inline double abs_max(int n, const double *x)
{ return std::fabs(x[index_abs_max(n, x)]); }

inline double abs_max(const std::vector<double> &x)
{ return std::fabs(x[index_abs_max(x)]); }

// saxpy (y=alpha*x+y) =======================================================


//inline void add_scaled(float alpha, const std::vector<float> &x, std::vector<float> &y)
//{ cblas_saxpy((int)x.size(), alpha, &x[0], 1, &y[0], 1); }

inline void add_scaled(float alpha, const std::vector<float> &x, std::vector<float> &y)
{
	size_t size = x.size();
	for(size_t i = 0; i < size; i++)
	{
		y[i] += alpha * x[i];
	}
}


};

#endif
